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Abstract 

Oscillatory behavior is a key property of many biological systems. The Small-Gain Theorem (SGT) 
for input/output monotone systems provides a sufficient condition for global asymptotic stability of an 
equilibrium and hence its violation is a necessary condition for the existence of periodic solutions. One 
advantage of the use of the monotone SGT technique is its robustness with respect to all perturba- 
tions that preserve monotonicity and stability properties of a very low-dimensional (in many interesting 
examples, just one-dimensional) model reduction. This robustness makes the technique useful in the 
analysis of molecular biological models in which there is large uncertainty regarding the values of kinetic 
and other parameters. However, verifying the conditions needed in order to apply the SGT is not al- 
ways easy. This paper provides an approach to the verification of the needed properties, and illustrates 
the approach through an application to a classical model of circadian oscillations, as a nontrivial "case 
study," and also provides a theorem in the converse direction of predicting oscillations when the SGT 
conditions fail. 
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1 Introduction 

Motivated by applications to cell signahng, our previous paper [I] introduced the class of monotone in- 
put/output systems, and provided a technique for the analysis of negative feedback loops around such 
systems. The main theorem gave a simple graphical test which may be interpreted as a monotone small 
gain theorem ("SGT" from now on) for establishing the global asymptotic stability of a unique equilibrium, 
a stability that persists even under arbitrary transmission delays in the feedback loop. Since that paper, 
various papers have followed-up on these ideas, see for example [261 IIZl HI 13 El Hffl 111 HI O [M] • The 
present paper, which was presented in preliminary form at the 2004 IEEE Conference on Decision and 
Control, has two purposes. 

The first purpose is to develop explicit conditions so as to make it easier to apply the SGT theorem, 
for a class of systems of biological significance, a subset of the class of tridiagonal systems with inputs 
and outputs. Tridiagonal systems (with no inputs and outputs) were introduced largely for the study of 
gene networks and population models, and many results are known for them, see for instance |3H I33j . 
Deep achievements of the theory include the generalization of the Poincare-Bendixson Theorem, from 
planar systems to tridiagonal systems of arbitrary dimension, due to Mallet-Paret and Smith [28] as well 
as a later generalization to include delays due to Mallet-Paret and Sell [27]. For our class of systems, we 
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provide in Theorem [T] sufficient conditions that guarantee the existence of characteristics ("nonUnear DC 
gain"), which is one of the ingredients needed in the SGT Theorem from [1]. 

Negative feedback is often associated with oscillations, and in that context one may alternatively 
view the failure of the SGT condition as providing a necessary condition for a system to exhibit periodic 
behaviors, and this is the way in which the SGT theorem has been often applied. 

The conditions given in Theorem [1] arose from our analysis of a classical model of circadian oscillations. 
The molecular biology underlying the circadian rhythm in Drosophila is currently the focus of a large 
amount of both experimental and theoretical work. The most classical model is that of Goldbeter, who 
proposed a simple model for circadian oscillations in Drosophila, see [BJ and the book [15j. The key to 
the Goldbeter model is the auto-inhibition of the transcription of the gene per. This inhibition is through 
a loop that involves translational and post-transcriptional modifications as well as nuclear translocation. 
Although, by now, several more realistic models are available, in particular incorporating other genes, see 
e.g. [241 125j . this simpler model exhibits many realistic features, such as a close to 24-hour period, and has 
been one of the main paradigms in the study of oscillations in gene networks. Thus, we use Goldbeter's 
original model as our "case study" to illustrate the mathematical techniques. 

The second purpose of this paper is to further explore the idea that, conversely, failure of the SGT 
conditions may lead to oscillations if there is a delay in the feedback loop. (As with the Classical Small Gain 
Theorem, of course the SGT is far from necessary for stability, unless phase is also considered.) As argued 
in [3], Section 3, and reviewed below, failure of the conditions often means that a "pseudo-oscillation" exists 
in the system (provided that delays in the feedback loop are large enough), in the rough sense that there 
are trajectories that "look" oscillatory if observed under very noisy conditions and for finite time-intervals. 
This begs the more interesting question of whether true periodic solutions exist. It turns out that some 
analogs of this converse result are known, for certain low-dimensional systems, see \29\ 122). In the context 
of failure of the SGT, Enciso recently provided a converse theorem for a class of cyclic systems, see 
The Goldbeter model is far from being cyclic, however. Theorem [2] in this paper proves the existence of 
oscillations for a class of monotone tridiagonal systems under delayed negative feedback, and the theorem 
is then illustrated with the Goldbeter circadian model. 

We first review the basic setup from [1]. 



2 I/O Monotone Systems, Characteristics, and Negative Feedback 

We consider an input /output system 

^ = f{x,u), y = h{x) , (1) 

in which states x{t) evolve on some subset X C M", and input and output values u{t) and y{t) belong 
to subsets U C R™ and y C RP respectively. The maps f : X x U ^ and h : X ^ Y are taken to 
be continuously differentiable. An input is a signal u : [0, oo) U which is locally essentially compact 
(meaning that images of restrictions to finite intervals are compact), and we write ip{t, xq,u) for the solution 
of the initial value problem dx/dt{t) = f{x{t),u{t)) with x(0) = xq, or just x{t) if xq and u are clear from 
the context, and y{t) = h{x{t)). Given three partial orders on X,U,Y (we use the same symbol -< for 
all three orders), a monotone input/output system ("MIOS"), with respect to these partial orders, is a 
system ([1]) which is forward-complete (for each input, solutions do not blow-up on finite time, so x{t) and 
y{t) are defined for all t > 0), /i is a monotone map (it preserves order) and: for all initial states xi,X2 for 
all inputs ui,U2, the following property holds: ii xi^X2 and ui<U2 (meaning that ui{t)^U2{t) for all t>0). 
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then ip(t, xi,u)^(p(t, X2,U2) for all t > 0. Here we consider partial orders induced by closed proper cones 
X C M^, in the sense that x^yiSy — x^K. The cones K are assumed to have a nonempty interior 
and are pointed, i.e. K f] —K = {0}. When there are no inputs nor outputs, the definition of monotone 
systems reduces to the classical one of monotone dynamical systems studied by Hirsch, Smith, and others 
[32], which have especially nice dynamics. Not only is chaotic or other irregular behavior ruled out, but, in 
fact, under additional technical conditions (strong monotonicity) almost all bounded trajectories converge 
to the set of steady states (Hirsch's generic convergence theorem [T9| [20]). 

The most interesting particular case is that in which K is an orthant cone in M", i.e. a set Sg of the 
form {x S M" | ffjXj > 0}, where = ±1 for each i. A useful test for monotonicity with respect to arbitrary 
orthant cones ("Kamke's condition" in the case of systems with no inputs and outputs) is as follows. Let 
us assume that all the partial derivatives ^{x,u) for i j, J^(x,u) for all and ^{x) for all i,j 
(subscripts indicate components) do not change sign, i.e., they are either always > or always < 0. We also 
assume that X is convex (much less is needed.) We then associate a directed graph G to the given MIOS, 
with n + m + p nodes, and edges labeled "+" or "— " (or ±1), whose labels are determined by the signs of 
the appropriate partial derivatives (ignoring diagonal elements of df /dx). One may define in an obvious 
manner undirected loops in G, and the parity of a loop is defined by multiplication of signs along the loop. 
(See e.g. [2| for more details.) A system is monotone with respect to some orthant cones in X, [/, Y if and 
only if there are no negative loops in G. In particular, if the cone is the main orthant (e = (1, . . . , 1)), the 
requirement is that all partial derivatives must be nonnegative, with the possible exception of the diagonal 
terms of the Jacobian of / with respect to x. A monotone system with respect to the main orthant is also 
called a cooperative system. This condition can be extended to non-orthant cones, see [30l [351 136 ] l37]. 

In order to define negative feedback ( "inhibitory feedback" in biology) interconnections we will say that 
a system is anti-monotone (with respect to given orders on input and output value spaces) if the conditions 
for monotonicity are satisfied, except that the output map reverses order: xi < X2 ^ h{x2) ^ h{xi). 

2.1 Characteristics 



A useful technical condition that simplifies statements (one may weaken the condition, see |26j) is that of 

existence of single-valued characteristics, which one may also think of as step-input steady-state responses 
or (nonlinear) DC gains. To define characteristics, we consider the effect of a constant input u{t) = mq, 
t > and study the dynamical system dx/dt = f{x,UQ). We say that a single- valued characteristic exists 
if for each uq there is a state K{uo) so that the system is globally attracted to K{uq), and in that case we 
define the characteristic k : U ^ y as the composition ho K. It is remarkable fact for monotone systems 
that (under weak assumptions on X, and boundedness of solutions) just knowing that a unique steady state 
K{uq) exists, for a given input value uq, already implies that K{uq) is in fact a globally asymptotically 
stable state for dx/dt = f{x,uo), see [23| [6]. 

2.2 Negative feedback 

Monotone systems with well-defined characteristics constitute useful building blocks for arbitrary systems, 
and they behave in many senses like one-dimensional systems. Cascades of such systems inherit the same 
properties (monotone, monostable response). Under negative feedback, one obtains non- monotone systems, 
but such feedback loops sometimes may be profitably analyzed using MIOS tools. 

We consider a feedback interconnection of a monotone and an anti- monotone input/output system: 




dxi 
~dF 



fl{xi,Ui), yi = hi{xi) 



(2) 
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-J^ = f2{x2,U2), y2 = h2{x2) (3) 

with characteristics denoted by "/c" and respectively. (We can also include the case when the second 
system is a static function y2it) = g{u2{t)).) As in [2j, we will require here that the inputs and outputs 
of both systems are scalar: mi=m2=pi=P2= 1; the general case [8] is similar but requires more notation 
and is harder to interpret graphically. The feedback interconnection of the systems ^ and ([3]) is obtained 
by letting n2=yi="y" and iti=y2="ii"; as depicted (assuming the usual real-number orders on inputs and 
outputs) in Figure [TJ 

characteristic of system 

— " — 

u y 





u 




characteristic of system 
\^ (anti-monotone) 




\, S 




N y 



Figure 1: Negative feedback configuration 

The main result from [1], which we'll refer to as the monotone SGT theorem, is as follows. We plot 
together k and g"^, as shown in Figure [21 and consider the following discrete dynamical system: 




Figure 2: characteristics 



= {g o k){u) 

on U . Then, provided that solutions of the closed-loop system are bounded, the result is that, if this 
iteration has a globally attractive fixed point u, as shown in Figure [2] through a "spiderweb" diagram, 
then the feedback system has a globally attracting steady state. (An equivalent condition, see Lemma 2.3 
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in [7], and [TIJ, is that the discrete system should have no nontrivial period-two orbits, i.e. the equation 
{g o k){u) = u has a unique solution.) 

It is not hard to prove, furthermore, that arbitrary delays may be allowed in the feedback loop. In other 
words, the feedback could be of the form u{t) = y{t—h), and such delays (even with h = h{t) time varying or 
even state-dependent, as long as t — h{t) ^ cxo as t ^ oo) do not destroy global stability of the closed loop. 
Moreover, it is also known [lOj that diffusion does not destroy global stability: a reaction-diffusion system, 
with Neumann boundary conditions, whose reaction can be modeled in the shown feedback configuration, 
has the property that all solutions converge to a (unique) uniform in space solution. 

2.3 Robustness 

It is important to point out that characteristics (called dose response curves, activity plots, steady-state 
expression of a gene in response to an external ligand, etc.) are frequently readily available from exper- 
imental data, especially in molecular biology and pharmacology, in contrast to the rare availability and 
high uncertainty regarding the precise form of the differential equations defining the dynamics and values 
for all parameters (kinetic constants, for example) appearing in the equations. MIOS analysis allows one 
to combine the numerical information provided by characteristics with the qualitative information given 
by "signed network topology" (Kamke condition) in order to predict global behavior. (See [34] for a longer 
discussion of this "qualitative-quantitative approach" to systems biology.) The conclusions from applying 
the monotone SGT are robust with respect to all perturbations that preserve monotonicity and stability 
properties of the 1-d iteration. 

Moreover, even if one would have a complete system specification, the 1-d iteration plays a role vaguely 
analogous to that of Nyquist plots in classical control design, where the use of a simple plot allows quick con- 
clusions that would harder to obtain, and be far less intuitive, when looking at the entire high-dimensional 
system. 

3 Existence of Characteristics 

The following result is useful when showing that characteristics exist for some systems of biological interest, 
including the protein part of the circadian model described later. The constant c represents the value of a 
constant control u{t) = c. 

Theorem 1 Consider a system of the following form: 



xq = c - ao{xo) + Po{xi) 



Xi = ai-i{xi-i) - Pi-i{xi) - ai{xi) + Pi{xi+i) 
i = 1, . . . ,n — 2 




+f]n-l{Xn) - 6{Xn-l) 
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evolving on M>q^, where c> is a constant. Assume that 9 and all the ai,(3i are differentiable functions 
[0, oo) [0, oo) with everywhere positive derivatives and vanishing at 0, 

and Oi, Pi, i = 0, . . . ,n — 2 are bounded, 

and 

an-i,f3n-i are unbounded. 

We use the notation 0{oo) to indicate lim^^oo ^(?'); CLnd similarly for the other bounded functions. Fur- 
thermore, suppose that the following conditions hold: 

ai_i((X)) + /3i(oo) < ai(oo) + /3j_i(oo) i = l,...,n-2 (4) 

^(oo) + /3i(oo) < aj(cx)), i = 0,...,n-2 (5) 
c < e{oo). (6) 
Then, there is a (unique) globally asymptotically stable equilibrium for the system. 

Observe that ([5]) (applied with i = 0) together with ([6]) imply that also: 

c + /3o(oo) < ao(oo) . (7) 

Proof. We start by noticing that solutions are defined for all t > 0. Indeed, consider any maximal solution 
x{t) = {xo{t),xi{t), . . .,Xn{t)). From 

^{xo + Xi + ...+Xn) = C-e{Xn^l) < C, (8) 

we conclude that there is an estimate Xi{t) < J2i Xi{t) < J2i Xi{0) + tc for each coordinate of x, and hence 
that there are no finite escape times. 

Moreover, we claim that x(-) is bounded. We first show that xq, . . . ,Xn-2 are bounded. For xq, it is 
enough to notice that xq < c — ao(xo) + f3o{oo), so that 

xo(t) > ao"^(c + ^o(oo)) =^ xo{t)<0 

so ([7|) shows that xq is bounded. Similarly, for Xj, i = 1, . . . , n — 2 we have that 

Xi < aj_i(cx)) - Pi-i{xi) - ai{xi) + Pi{oo) 

so dH) provides boundedness of these coordinates as well. 
Next we show boundedness of Xn—i and Xn- 

Since the system is a strongly monotone tridiagonal system, we know (see [31j . Corollary 1), that x„(i) 
is eventually monotone. That is, for some T > 0, either 

Xnit) > Vt > T (9) 

or 

Xn{t) < Vt > T. (10) 

Hence, x„(t) admits a limit, either finite or infinite. 
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Assume first that Xn is unbounded, which means that oo because of eventual monotonicity. 

Then, case (fTOjl cannot hold, so ([9]) holds. Therefore, 



a„_i(x„_i(t)) - /3„_i(x„(t)) = Xn > 

for all t >T, which implies that 

Xn-i{t) > a~li(/3„_i(x„(t))) ^ cx) 
as well. Looking again at ([8]), and using that c — 6{oo) < (property dS])), we conclude that 

^ (Xo + Xi + . . . + Xn-l + Xn) (t) < 

for all t sufficiently large. Thus xq + xi + . . . + Xn~i + x„ is bounded (and nonnegative) , and this implies 
that Xn-i is bounded, a contradiction since we showed that — > oo. So Xn is bounded. 

Next, notice that Xn-i < an_2(x„_2) + Pn-i{xn) — Oin-i{xn-i) ■ The two positive terms are bounded, 
because both x„_2 and bounded. Thus, 

Xn-l < v{t) - an-l(x„_i) , 

where < v{t) < k for some constant k. Thus x„_i(t) < whenever x„_i(t) > a~l;^(/c), and this proves 
that Xn-l is bounded, as claimed. 

Once that boundedness has been established, if we also show that there is a unique equilibrium, then the 
theory of strongly monotone tridiagonal systems I32j) (or [221 E] ^^r more general monotone systems 
results) will ensure global asymptotic stability of the equilibrium. So we show that equilibria exist and are 
unique. 

Let us write /j(x) for the right-hand sides of the equations, so that Xj = fi{x) for each i. We need to 
show that there is a unique nonnegative solution x = (xq, . . . , Xn) of 

fo{x) = . . . = fn{x) = 0. 

Equivalently, we can write the equations like this: 

fn{x) = (11) 
n{x) + . . . + fnix) = (12) 

/o(x) + /i(x)... + /„(x) = (13) 

Since /o(x) + /i(x) . . . + /n(x) = c - 6{ Xn— i), (jlSp has the unique solution Xn—i — Xn—i — ^ ^(c) which is 
well-defined because property ([6|) says that c < 9{oo). 

Next, we consider equation pT|) . This equation has the unique solution: 

Xn — Xn — Pn—li'^n—liXn—l')) 

which is well-defined because Pn-i is a bijection. 
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Pick i E {1, . . . , n — 1} and suppose that we have uniquely determined xj = Xj for each j > i. We will 
show that Xi-i is also uniquely defined. Equation ([12]) is: 



ai-i{xi^i) - Pi^i{xi) - 9{xn-i) = 

and has the unique solution 

Xi-i = Xi-i = a~\{(3i-i{xi) + 6{xn-i)) 

which is well-defined because property ([5]) says that 9{oo) + (3i-i{oo) < ai^i{oo) for each i = 1, . . . ,n — 1. 
By induction on j = n — 1, . . . , 1, we have completed the uniqueness proof. I 



4 The Goldbeter Circadian Model 



The original Goldbeter model of Drosophila circadian rhythms is schematically shown in Figure [3l The 



M-h^ 



V2 



V4 



ki 



k2 



Figure 3: Goldbeter's Model 



assumption is that PER protein is synthesized at a rate proportional to its mRNA concentration. Two 
phosphorylation sites are available, and constitutive phosphorylation and dephosphorylation occur with 
saturation dynamics, at maximum rates Vi and with Michaelis constants K^. Doubly phosphorylated PER 
is degraded, also described by saturation dynamics (with parameters v^, kd), and it is translocated to the 
nucleus, with rate constant ki. Nuclear PER inhibits transcription of the per gene, with a Hill-type reaction 
of cooperativity degree n and threshold constant Kj. The resulting mRNA is produced, and translocated 
to the cytoplasm, at a rate determined by a constant Vs- Additionally, there is saturated degradation of 
mRNA (constants and km)- 

Corresponding to these assumptions, and assuming a well-mixed system, one obtains an ODE system 
for concentrations are as follows: 



M 



+ km + M 

k.M-^^+ ^^^^ 



ViPo V2P1 VsPi ^ V4P2 ^^^^ 



Ki+Po K2 + P1 K3 + P1 K4 + P2 
p V3P1 V4P2 , p ^ , p VdP2 

P2 = — — — - kiP2 + K2 "at 



7^3+ A K4 + P2 kd + P2 

Pn = klP2-k2PN 

where the subscript i = 0,1,2 in the concentration Pj indicates the degree of phosphorylation of PER 
protein, Pjy is used to indicate the concentration of PER in the nucleus, and M indicates the concentration 
of per mRNA. 
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Parameter 


Value 


Parameter 


Value 


k2 


1.3 


ki 


1.9 


Vi 


3.2 


V2 


1.58 


Vz 


5 


Vi 


2.5 


Vs 


0.76 


km 


0.5 


kg 


0.38 


Vd 


0.95 


kd 


0.2 


n 


4 


Ki 


2 


K2 


2 


K3 


2 


K4 


2 


Kj 


1 


Vm 


0.65 



Table 1: Parameter Values 

The parameters (in suitable units fiM or h~^) used by Goldbeter are given in Table [H With these 
parameters, there are limit cycle oscillations. If we take Vg as a bifurcation parameter, a Hopf bifurcation 
occurs at Vg ~ 0.638. 

As an illustration of the SGT, we will show now that the theorem applies when Vg = 0.4. This means 
that not only will stability of an equilibrium hold globally in that case, but this stability will persist even 
if one introduces delays to model the transcription or translation processes. (Without loss of generality, 
we may lump these delays into one delay, say in the term appearing in the equation for M.) On the 
other hand, we'll see later that the SGT discrete iteration does not converge, and in fact has a period-two 
oscillation, when Vg = 0.5. This suggests that periodic orbits exist in that case, at least if sufficiently large 
delays are present, and we analyze the existence of such oscillations. 

For the theoretical developments, we assume from now on that 

Vg < 0.54 (15) 

and the remaining parameters will be constrained below, in such a manner that those in the previously 
given table will satisfy all the constraints. 

4.1 Breaking-up the Circadian System and Applying the SGT 

We choose to view the system as the feedback interconnection of two subsystems, one for M and the other 
one for P, see Figure HI 



Ul 


M 


Vi 














p 




y2 





Figure 4: Systems in feedback 
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mRNA Subsystem 



The mRNA (M) subsystem is described by the scalar differential equation 



with input ui and output yi = kgM. 

As state-space, we will pick a compact interval Xi = [0,M], where 



<M<^ (16) 



I'm ~ "^s k^ 

and we assume that Vg < v^- The order on Xi is taken to be the usual order from R. 
Note that the first inequality implies that 



and therefore 



VmM 

V, < ^ (17) 



VsK? VmM 
SI m ^ Q 



Kf + u'l km + M 
for all ui > 0, so that indeed Xi is forward-invariant for the dynamics. 

With the parameters shown in the table given earlier (except for Vg, which is picked as in psp ). 

M = 2.45 

satisfies all the constraints. 

As input space for the mRNA system, we pick Ui = M>o, and as output space Yi = [0,^^). Note that 
yi = kgM < kgM < v^, by pUj) . so the output belongs to Yi. We view Ui as having the reverse of the 
usual order, and Yi are is given the usual order from M. 

The mRNA system is monotone, because it is internally monotone {df/du < 0, as required by the 
reverse order on Ui) and the output map is monotone as well. 

Existence of characteristics is immediate from the fact that M > for M < k{ui) and M < for 
M > k{ui), where, for each constant input ui, 

M Vg Kf km 
ui) - 



VmKf + VmUi^-VgKf 

(which is an element of Xi). 

Note that all solutions of the differential equations which describe the M-system, even those that do 
not start in Xi, enter Xi in finite time (because M(t) < whenever M(t) > M, for any input mi(-)). The 
restriction to the state space Xi (instead of using all of M>o) is done for convenience, so that one can view 
the output of the M system as in input to the P-subsystem. (Desirable properties of the P-subsystem 
depend on the restriction imposed on U2-) Given any trajectory, its asymptotic behavior is independent 
of the behavior in an initial finite time interval, so this does not change the conclusions to be drawn. 
(Note that solutions are defined for all times -no finite explosion times- because the right-hand sides of 
the equations have linear growth.) 
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Protein Subsystem 



The second (P) subsystem is four-dimensional: 

Pq = U2- — — —=- + 



Ki + Po K2 + P1 
ViPo V2P1 V^Pi , F4P2 



1 



Ki+Po K2 + P1 K3 + P1 K4 + P2 

"2 = T7 — T7 — hP2 + K2 "Af 



i^3+Pl K^ + P2 kd + P2 

Pn = hP2-k2PN 



with input U2 and output 1/2 = Pn- 

For the P subsystem, the state space is R>g with the main orthant order, and the input space is 
U2 = Yi and the output space is Y2 = Ui (with the orders specified earher). Internal monotonicity of 
the P subsystem is clear from the fact that |^ > for all i / j (cooperativity) . In fact, because these 
inequalities are strict and the Jacobian matrix is tridiagonal and irreducible at every point, this is an 
example of a strongly monotone tridiagonal system ( \31\ I32j). The system is anti-monotone because the 
identity mapping reverses order (recall that Y2 = Ui has the reverse order, by definition). 

We obtain the following result as a corollary of Theorem [U applied with n = 3, 0(r) = .'"'^ , ao(^) = 
j^^^j. , etc. It says that, for the parameters in Table [H as well as for a larger set of parameters, the system 
has a well-defined characteristic, which we will denote by g. (It is possible to give an explicit formula for 
g~^ , in this example.) 

Proposition 4.1 Suppose that the following conditions hold: 

• Vd + V2<Vi 

• Vi+V4<V2 + V3 

• < c < Vd 

• VA + Vd<V3 

and that all constants are positive and the input U2{t) = c. Then the P-system has a unique globally 
asymptotically stable equilibrium. 



5 Closing the Loop 

Solutions of the closed-loop system, i.e., of the original system (fH|) . are bounded under the above assump- 
tions. To see this, we argue as follows. Take any solution of the closed loop system. As we pointed out 
earlier, there are no finite time explosions, and also the M coordinate will converge to the set Xi = [0, M]. 

This means that the subsystem corresponding to the P-coordinates will be forced by an input U2 such 
that u{t) £ [0, kgM] for all t > tQ, for some tQ. Now, for constant inputs in [0, v^), which contains [0, kgM], 
we have proved that a characteristic k exists for the open-loop system corresponding to these coordinates. 
Therefore, by monotonicity, the trajectory components y{t) = (Po(t), Pi(i), P2(i), PAr(t)) will lie in the 
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main orthant order rectangle [yo(0> fo'^ each t > 0, where yo is the solution with constant input 
U2 = and yo{to) = y{to) and where y^ is the solution with constant input U2 = kgM, and yi(io) = y(*o)- 
Since yo and yi converge to [k{0), k{ksM)], the omega-limit set of y is included in [k{0), k{ksM)], and 
therefore the P components are bounded as well. 

Now we are ready to apply the main theorem in p|. In order to do this, we first need to plot the 
characteristics. See Figure [5] for the plots of g and k~^ (dashed and dotted curves) and the a typical 
"spiderweb diagram" (solid lines), when we pick the parameter Vs = 0.4. It is evident that there is global 



Figure 5: Stability of spiderweb {vg = 0.4) 

convergence of the discrete iteration. Hence no oscillations can arise, even under arbitrary delays in the 
feedback from P/v to M, and in fact that all solutions converge to a unique equilibrium. 

On the other hand, for a larger value of Vg, such as Vs = 0.5, the discrete iteration conditions are 
violated; see Figure [6] for the "spiderweb diagram" that shows divergence of the discrete iteration. Thus, 



Figure 6: Instability of spiderweb {vg = 0.5) 
one may expect periodic orbits in this case. We next prove a result that shows that indeed that happens. 

6 Periodic Behavior when SGT Conditions Fail 

One may conjecture that there is a connection between periodic behaviors of the original system, at least 
under delayed feedback, and of the associated discrete iteration. We first present an informal discussion 
and then give a precise result. 

For simplicity, let us suppose that k already denotes the composition of the characteristics g and k. 
The input values u with k{k{u)) = u which do not arise from the unique fixed point of k{u) = u are 
period- two orbits of the iteration = k{u). Now suppose that we consider the delay differential system 
x{t) = f{x(t), h{x(t — r))), where the delay r > is very large. We take the initial condition x{t) = xq, 
t £ [— r, 0], where xq is picked in such a manner that h{xo) = uq, and uq ^ u\ are two elements of \J such 
that kiuo) = u\ and k{u\) = uq- If the input to the open-loop system x = f{x,u) is u{t) = uq, then the 
definition of characteristic says that the solution x{t) approaches xi, where h{xi) = ui, Thus, if the delay 
length r is large enough, the solution of the closed-loop system will be close to the constant value xi for 
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t r. Repeating this procedure, one can show the existence of a Hghtly damped "oscillation" between the 
values xq and xi, in the sense of a trajectory that comes close to these values as many times as desired 
(a larger r being in principle required in order to come closer and more often). In applications in which 
measurements have poor resolution and time duration, it may well be impossible to practically determine 
the difference between such pseudo-oscillations and true oscillations. 

It is an open question to prove the existence of true periodic orbits, for large enough delays, when the 
small-gain condition fails. The problem is closely related to questions of singular perturbations for delay 
systems, by time-re-parametrization. We illustrate this relation by considering the scalar case, and with 
y = X. The system x = f{x,x{t — r)), has periodic orbits for large enough r if and only if the system 
£x{t) = f{x{t),x{t — 1)) has periodic orbits for small enough e > 0. For e = 0, we have the algebraic 
equation f{x,u) = that defines the characteristic x = k{u). Thus one would want to know that periodic 
orbits of the iteration = k(u), seen as the degenerate case e = 0, survive for small e > 0. A variant of 
this statement is known in dimension one from work of Nussbaum and Mallet-Paret (|29j). which shows 
the existence of a continuum of periodic orbits which arise in a Hopf bifurcation and persist for < e <C 1; 
see also the more recent work [22]. (We thank Hal Smith for this observation.) 

We now show that, at least, for a class of systems which is of some general interest in biology, and 
which contains the circadian model, oscillations can be proved to exist if delays are large enough and the 
SGT fails locally (exponential instability of the discrete iteration). 

6.1 Predicting Periodic Orbits when the Condition Fails 

In this section we prove the following theorem, which applies immediately to the complete circadian 
model (fni) . 

Theorem 2 Consider a tridiagonal system x = f{x,u) with scalar input u: 



and scalar output y = Xn- The functions fi are twice continuously differentiable, and (cooperativity) all the 
off-diagonal Jacobian entries are positive. Suppose that there is a unique pair {xq,uq) G x M such that 
f{xo,uo) = 0, and consider the linearized system {A,b,c), where b = (1, 0, . . . , 0, 0)', c = (0, 0, . . . , 0, 1), 
and A = D^f, the Jacobian of the vector field f{x, uq) evaluated at xq. Assume that A is nonsingular, and 
let g = c{—A)~^b be the DC gain of the linearized system. Pick any positive number k such that kg > 1, 
and consider a delayed feedback u{t) = —ky(t — h) with h > 0. Then, for some h > 0, the system U8\) under 
the feedback u{t) = —ky{t — Hq) admits a periodic solution, and, moreover, the omega-limit set of every 
bounded solution is either a periodic orbit, the origin, or a nontrivial homoclinic orbit with lim^^-i-oo = a^o- 

Note that the uniqueness result for closed-loop equilibria will always hold in our case, and the DC gain 
property kg > 1 corresponds to a locally unstable discrete iteration. The matrix A is Hurwitz when we have 
hyperbolicity and parameters as considered earlier (existence of characteristics). The conclusion is that, 
for a suitable delay length /iq, there is at least one periodic orbit, and, moreover, bounded solutions not 
converging to zero exhibit oscillatory behavior (with periods possibly increasing to infinity, if the omega- 
limit set is a homoclinic orbit). (We conjecture that moreover, for the circadian example, in fact almost all 



X2 = 



Xl = 



fl{xi,X2,u) 
f2{xi,X2,Xs) 



(18) 
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solutions converge to a periodic orbit. Proving this would require establishing that no homoclinic orbits 
exist for our systems, just as shown, when no delays present, in |28|.) 

Before proving Theorem [2l we show the following simple lemma about linear systems. 



Lemma 6.1 Consider a linear n-dimensional single-input single-output system {A, b, c), with b = (1, 0, . . . , 0, 
and c = (0, 0, . . . , 0, 1), and suppose that j4 is a linear tridiagonal matrix 



A = 





h2 










d2 


^3 








as 


^3 


64 


Vo 












with Ojftj > for all i (in particular, this holds if all off-diagonal elements are positive). Then, the transfer 
function W{s) = c{sl — A)~^b has no zeroes and has distinct real poles; more specifically, W{s) = 
where Po = 02 ... (in and q{s) = (s— ai) . . . (s— a„) for distinct real numbers ai, . . . , a„. Moreover, there are 
two real- valued functions /i : C — > M and v : C ^ R>o so that the logarithmic derivative Q{s) = q'{s)/q{s) 
satisfies Q{s) = iJ,{s) — ii^(s)Ims for every s that is not a root of q. 



Proof. The fact that A has n distinct real eigenvalues is a classical one in linear algebra; we include a short 
proof to make the paper more self-contained. Pick any positive number ai and define, inductively, 

for i = 2, . . . ,n. Let S = diag (di, . . . , cr„). Then B = S^^AS is a tridiagonal symmetric matrix: 

/di C2 ... \ 
C2 d2 C3 ... 
^ ^ C3 ds C4 ... 

V ... Cn dn) 

where Cj = ei\/aibi and = signoj = sign 6 j G {— 1,+1}. Therefore B, and hence also A, has all its 
eigenvalues real. Moreover, there is a basis {vi, . . . ,Vn} consisting of orthogonal eigenvectors of B, and 
so A admits the linearly independent eigenvectors Svi. Moreover, all eigenvalues of B (and so of A) are 
distinct. (Pick any A and consider C .= B — XI. The first n — 1 rows of C look just like those of B, with 
di := di — X. The n — 1 x n matrix consisting of these rows has rank n — 1 (just consider its last n — 1 
columns, a nonsingular matrix), so it follows that C has rank> n — 1. Therefore, the kernel of C has 
dimension at most one.) We conclude that A has n distinct real eigenvalues and hence its characteristic 
polynomial has the form q{s) = (s — ai) . . . (s — a„). 

By Cramer's rule, {si — A)~^ = (l/q(s))cof(s/ — A), where "cof" indicates matrix of cofactors. Thus 
W{s) = Po/q{s), where pq is the (n, 1) entry of cof(s/ — A), i.e. (— l)""*"-^ times the determinant of the 

matrix {si — A)-^ ^ obtained by deleting from si — A the first row and last column. The matrix {si — A)^ ^ 
is upper triangular, and its determinant is (—02) . . . {—an) = (— l)"~^a2 • • • ^n- Therefore po = 02 • • • an, as 
claimed. 
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Finally, consider Q{s) = Ylk=i s+aj, • Write s = a + ib, so that 

1 1 (a + Ok) — ib 



s + ak (a + Ok) + ib pk 
where pk = (a + a^)^ + 6^, and therefore 



s) — iv{s)b 



as desired. 



We now continue the proof of Theorem [21 by first studying the closed- loop linearized system x{t) 
Ax{t) + bkcx{t — h). The closed-loop transfer function 

W 



1 + ke-^'>W{s) 

corresponding to a negative feedback loop with delay h and gain k, simplifies to: 

Wh = ^, F{s,h)=qis)+pe~^^, 

where p = pok. 

In order to prove that there are oscillatory solutions for some h = Hq, we proceed as follows. We 
will use the weak form of the Hopf bifurcation theorem ( "weak" in that no assertions are made regarding 
super or subcriticality of the bifurcation) as given in Theorem 11.1.1 in [T8]. The theorem guarantees that 
oscillatory solutions will exist, for the nonlinear system, and for some value of the delay h arbitrarily close 
to a given /iq > 0, provided that the following two properties hold for ho: 

(HI) There is some ljq ^ such that F{iujQ,ho) = 0, w = iujQ is a simple root of F{uj,ho) = 0, and 
(nonresonance) F{miuJ,hQ) ^ for all integers m > 1; 

and letting A(/i) be a function such that F{X{h), h) = for all h near ho and A(/io) = i^o (such a function 
always exists): 

(H2) ReA'(/io) / 0. 

In order to prove these properties, we proceed analogously to what is done for cyclic systems in [9]. 
(Cyclic systems are the special case in which dfi/dxi+i = for each i = l,...,n — 1, which is not the case 
in our circadian system.) 

We first show that F{iu>o,ho) = for some ho > and uio > 0. Since q{s)F{s,h) = 1 -|- {p/q{s))e~^^ 
and q{iio) ^ for all real numbers uj (because q has only real roots, and A is nonsingular, so also q{0) ^ 0), 
it is enough to find an /iq > and > such that /(wq) = —e^'^°'^° = e^i^o^^o--^) ^ where /(w) = p/q{iu)). 
Since / is a continuous function on [0,oo), /(O) = Pok/q{0) = W{0)k = gk > 1 hy assumption, and 
lim^^^oo /(^) = 0, there is some > such that |/(w)| = 1, so that /(w) = e**^ for some f which we may 
take in the interval (0, 2tt]. It thus suffices to pick /iq = + vr /wq , so that /iqWo — tt = (p. 

Fix any such ho- Since for retarded delay equations there are at most a finite number of roots on 
any vertical line, we can pick ujo with largest possible magnitude, so that F{miuJo, ho) 7^ for all integers 
m > 1. To prove that (HI) and (H2) hold for these ho and wq, we first prove that {dF/ds){iuj, h) is nonzero 
at each point iu in the imaginary axis, and all > (so, in particular, at (za;o,/io))- By the implicit 
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function theorem, this wih imply that is a simple root, as needed for (HI). Since F{X{h),h)) = in a 
neighborhood of ho, we also have that 

dF BF 
A'W = -^{Ks),h)/—(Ks),h) 

spe "■^ s 



q'{s) - hpe-^" _ /j ■ 

p^ — hs 

At o; = u;o and h = /iq, we have that F{u)^h) = 0, so ^(5) = — pe~^^, and therefore, denoting Q{s) 
q'{s)/q{s): 



Q{iuJo) + ho ' 



It follows that, for ujo > 0: 



ReX'iho) = Im ^ 



Q{iLo) + ho 
= -Im [Q{iuJo) + ho] = -lmQ{iuJo) ■ 

The proof that periodic orbits exist, for each h near enough ho, will then be completed if we show that 
lmQ{iu;) 7^ for all nonzero real numbers uj. But, indeed. Lemma l6.1l savs that Q{iuj) = fj,{iuj) — iu{iuj)ijo, 
where > 0, so linQ{iuj) = —v{iijj)uj 7^ 0. 

To conclude the proof, we note that the conclusion about global behavior follows from the Poincare- 
Bendixson for delay-differential tridiagonal systems due to Mallet-Paret and Sell [27]. I 

Note that, since ReA'(/io) 7^ 0, if /iq is near enough ho, then the system (fTSj) under negative feedback 
u = —ky{t — /iq) admits a pair of complex conjugate eigenvalues a + iuj for its linearization, with a > 0. 
Thus, its equilibrium is exponentially unstable, and therefore every bounded solution not starting from the 
center-stable manifold will in fact converge to either a homoclinic orbit involving the origin or a periodic 
orbit. 



6.2 Examples 



As a first example, we take the system with the parameters that we have considered, and Vg = 0.5. We have 
seen that the spider-web diagram suggests oscillatory behavior when delays are present in the feedback 
loop. We first compute the equilibrium of the closed- loop system (with no delay), which is approximately: 



M 1.47, Pi = 0.42, P2 = 0.29, Pq = 0.71, Pj 



N 



0.42. 



We now consider the system with variables M, Pi, P2, P3, P/v in which the feedback term VgKf /{Kf + Pg^) 
is replaced by an input u. Let A be the Jacobian of this open-loop dynamics evaluated at the positive 
equilibrium given above. Then si — A ^ 



0.08 + s 
-0.38 








0.87 + s 
-0.87 







-0.54 
2.24 + s 
-1.70 







-0.96 
3.66 + s 
-1.9 







-1.3 
1.3 + s 
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and hence the transfer function W{s) = c{sl — A) ^b, where b = (1, 0, 0, 0, 0)' and c = (0, 0, 0, 0, 1), is: 



Wis) = 




where po « 1.075 and q{s) 

(0.084 + s) (s^ + 8.08 + 17.61 + 10.98 s + 1.56) . 

The DC gain of the system is g = W{0) « 8.26 (which is positive, as it should be since the open loop system 
is monotone and has a well-defined steady-state characteristic), and k = —d{vsKf/(Kf + P^)/dP2, ~ 
0.14 when evaluated at the computed equilibrium. Thus /(O) = gk k, 1.138 > 1, as required. Indeed, 
Im. Q{iijj) 

- (2.88 + 133.26a;2 + 408.07 o;"^ + 120.12tj6 + S.Ow^) uj 
(0.007 + (2.42 + 65.68lj2 + 135.81^^ _^ 3o.02cj6 + ^^8) 

and hence ImQ(zcj) ^ for all oj ^ 

We show in Figure [7] one simulation, with h = 100, showing a periodic limit cycle. The delay length 




1000 1100 1200 1300 1400 1500 1600 1700 1S00 1900 2000 



Figure 7: Oscillations seen in simulations [vg = 0.5, delay = lOOhr, initial conditions all at 0.2), using 
MATLAB's dde23 package 

needed for oscillations when Vg = 0.5 is biologically unrealistic, so we also show simulations for V2 = 0.6, a 
value for which no oscillations occur without delays, but for which oscillations (with a period of about 27 
hours) occur when the delay length is about 1 hour, see Figure El 

7 A Counterexample 

We now provide a (non- monotone) system as well as a feedback law u = g{y) so that: 

• the system has a well-defined and increasing characteristic k, 

• the discrete iteration ti+ = g{k{u)) converges globally, and solutions of the closed-loop system are 
bounded, 

yet a stable limit-cycle oscillation exists in the closed-loop system. This establishes, by means of a simple 
counterexample, that monotonicity of the open-loop system is an essential assumption in our theorem. 
Thus, robustness is only guaranteed with respect to uncertainty that preserves monotonicity of the system. 
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Figure 8: Oscillations seen in simulations {vg = 0.6, delay = Ihr, initial conditions all at 0.2), using 
MATLAB's dde23 package 



The idea underlying the construction is very simple. The open-loop system is linear, and has the 
following transfer function: 

Wis) = ^^±i 

^ ' s2 + (0.25)s + 1 

Since the DC gain of this system is PF(0) = 1, and the system is stable, there is a well-defined and increasing 
characteristic k{u) = u. However a negative feedback gain of 1/2 destabilizes the system, even though the 
discrete iteration u'^ = (— l/2)n is globally convergent. (The Hoo gain of the system is, of course, larger 
than 1, and therefore the standard small-gain theorem does not apply.) In state-space terms, we use this 
system: 

Xl = {—l/4:)xi—X2 + 2u 

X2 = Xl 

y = (l/2)(x2-xi). 

Note that, for each constant input u = uq, the solution of the system converges to (0, no/2), and therefore 
the output converges to uq, so indeed the characteristic k is the identity. 

We only need to modify the feedback law in order to make solutions of the closed-loop globally bounded. 
For the feedback law we pick g{x) = — 0.5sat(y), where sat(-) := sign(-) min{l, | • |} is a saturation function. 
The only equilibrium of the closed-loop system is at (0,0). 

The discrete iteration is 

u+ = -(l/2)sat(u) . 

With an arbitrary initial condition uq, we have that ui = — (l/2)sat(no), so that < 1/2. Thus 
Uk = {—l/2)uk-i for all k >2, and indeed — > so global convergence of the iteration holds. 

However, global convergence to equilibrium fails for the closed-loop system, and in fact there is a 
periodic solution. Indeed, note that trajectories of the closed loop system are bounded, because they can 
be viewed as solutions of a stable linear system forced by a bounded input. Moreover, since the equilibrium 
is a repelling point, it follows by the Poincare-Bendixson Theorem that a periodic orbit exists. Figure [9] is 
a simulation showing a limit cycle. 
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Figure 9: Limit cycle in counterexample 
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